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Abstract. Large linear systems with sparse, non-symmetric matrices arise in the modeling of 
Markov chains or in the discretization of convection-diffusion problems. Due to their potential 
to solve sparse linear systems with an effort that is linear in the number of unknowns, algebraic 
multigrid (AMG) methods are of fundamental interest for such systems. For symmetric positive 
definite matrices, fundamental theoretical convergence results are established, and efficient AMG 
solvers have been developed. In contrast, for non-symmetric matrices, theoretical convergence 
results have been provided only recently. A property that is sufficient for convergence is that 
the matrix be an M-matrix. In this paper, we present how the simulation of incompressible fluid 
flows with particle methods leads to large linear systems with sparse, non-symmetric matrices. 
In each time step, the Poisson equation is approximated by meshfree finite differences. While 
traditional least squares approaches do not guarantee an M-matrix structure, an approach based 
on linear optimization yields optimally sparse M-matrices. For both types of discretization 
approaches, we investigate the performance of a classical AMG method, as well as an AMLI 
type method. While in the considered test problems, the M-matrix structure turns out not to 
be necessary for the convergence of AMG, problems can occur when it is violated. In addition, 
the matrices obtained by the linear optimization approach result in fast solution times due to 
their optimal sparsity. 



1. Introduction 

The efficient and robust numerical approximation of the Poisson equation is a key problem in many 
applications. In some special cases (such as spectral methods that use the fast Fourier transform), 
the Poisson equation can be solved in a single step. However, in most scenarios, the solution 
process is done in two consecutive steps: first, a discretization transforms the Poisson equation 
into a finite linear system; second, the linear system is solved. While in most investigations on 
multigrid methods the discretization is assumed as canonically given, in this paper we present 
an application in which many possible discretization strategies exist, and the interplay between 
discretization and performance of multigrid solvers is of interest. 

We consider Poisson problems that arise in Lagrangian particle methods for the simulation of 
incompressible fluid flows. Particle methods solve the governing fluid equations on a set of nodes 
that move with the flow. Thus the discretization is done in the more natural Lagrangian frame of 
reference. As a consequence, convective terms are solved exactly. One of the earliest Lagrangian 
particle methods is smoothed particle hydrodynamics (SPH) |29[ 118] . which was first applied to 
astrophysical fluid dynamics, and later expanded to a wider class of hydrodynamic equations 
[32|, I33j . Generalization of the SPH approach, that use moving least squares approximations for 
derivatives, have been presented by Dilts [12] and Kuhnert [22]. Both approaches have been 
successfully applied to problems with complex, time-dependent geometries [13] and free surface 
flows [24| . Further examples are shown in [44J . The problems considered in this paper are mostly 
formulated in the terminology of those latter two approaches. However, many aspects transfer 
to other types of particle methods. A fundamental difficulty with Lagrangian particle methods 
is that one gives away the control over the positions of the computational nodes. While one can 
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(and in fact should) remove particles where not desired, and insert new particles where required, 
the points are in general completely unstructured and without connections between each other. 

A classical solution strategy for incompressible flows is the Chorin projection method [7], which 
has been successfully applied to Lagrangian particle methods miTO]. In each time step, the flow 
is first evolved neglecting pressure, which yields a velocity field that is not incompressible. Then 
the velocity field is modified by a gradient field that is chosen so that incompressibility is restored. 
This correction step requires the solution of a Poisson equation. As noted above, the arrangement 
of the particles is in a general, unstructured cloud of points. And the Poisson equation has to be 
approximated on this geometry. 

One possibility to approximate the Poisson equation is using a finite element approach. For 
instance, one can generate a triangulation [46] (or another type of mesh) between the points, and 
apply a traditional finite element approach [47j . However, mesh generation is an expensive task, 
especially in 3d, and a good mesh generator may not always be available. In the context of particle 
methods, one rarely is willing to make such an investment. An alternative is the use of meshfree 
finite element approaches, such as the element-free Galerkin method (EFGM) [5], or the partition 
of unity finite element method (PUFEM) [31 [3|. Similar approaches [19] have been applied in 
a multilevel context [20]. A key advantage of all these approaches is that they are based on 
inner products between (meshfree) basis functions, hence they lead to symmetric linear systems 
to be solved. A fundamental drawback is that the computation of the inner products requires 
quadratures, which in many cases can only be achieved by introducing a regular background grid. 
This step is not only technically inconvenient, it is also costly. Another difficulty with finite 
element approaches is the treatment of general boundary conditions, which arise in multi-phase 
flows in complex, time-dependent domains. 

In contrast, meshfree finite difference methods do not require basis functions or quadrature, and 
boundary conditions can be implemented directly. At each point, the approximation of the Laplace 
operator is based on a local Taylor expansion. A first approach by Perrone and Kao [36| was later 
extended by Liszka, Orkisz, et al. [281 111] to the generalized finite difference method (GFDM). A 
detailed presentation of common particle methods and meshfree approaches is provided in j43j . In 
this paper, we outline the fundamental conditions for a consistent approximation of the Laplace 
operator in Sect. [21 

A key drawback of meshfree finite difference approaches is that the resulting linear systems are in 
general non-symmetric. As we explain in Sect. [3l this non-symmetry cannot be overcome by any 
simple means. In the same section, we outline theoretical results on the convergence of algebraic 
multigrid approaches for non-symmetric matrices. A property sufficient for convergence is an 
M-matrix structure. We present fundamental properties of M-matrices and their relevance for 
multilevel methods in Sect.[4l 

In addition to the non-symmetry of the meshfree finite difference matrices classical least squares 
approaches, as outlined in Sect. [HI involve two other difficulties. First, the arising matrices are 
about six times as dense as finite difference matrices on rectangular grids. Second, an M-matrix 
structure is not guaranteed, and is — unless the point geometry is sufficiently nice — in general vi- 
olated. Unlike the non-symmetry, these latter two problems can be overcome, by a new approach 
that is based on linear sign-constrained optimization. This approach is presented in Sect. [3 and 
conditions for its successful application are provided in Sect. [51 In Sect. [71 theoretical predic- 
tions about the computational cost to generate the Poisson matrices with least squares vs. linear 
optimization approaches are given. In Sect. [H we numerically investigate the performance of 
various multigrid approaches for the matrices constructed by the different types of minimization 
approaches. Conclusions and an outlook are given in Sect.[9l 
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2. Meshfree Finite Differences for the Poisson Equation 

Consider the Poisson equation to be solved inside a domain il C K.'^ 

'Aw — f in Q 



9 onTo (1) 



u ■ 



^ = h on Fat 
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where F^i U Fat = dVl. We consider the case of ([T]) possessing a unique smooth solution. 

Fet a point cloud X — {xi, . . . , x„} C r2 be given, which consists of interior points Xi <Z U, and 
boundary points Xb C d^l. The point cloud is meshfree, i.e. no information about connection of 
points is provided. Meshfree finite difference approaches convert problem ([1]) into a linear system 

A-i=f, (2) 

where the vector u contains approximations to the values u{xi). The i-th row of the matrix A 
consists of the stencil corresponding to the point Xi . 

In the following, we derive conditions on finite difference stencils to yield an approximation to 
the Laplace operator. Consider a function u G C'^{il C At each interior point xq, we 

wish to approximate AM(a;o) using the function values of a finite number of points in a circular 
neighborhood {xo,xi, . . . ,Xm) G B{xo,r), where B{xQ,r) = {x & fl : \\x ~ xo\\2 < r}. Note that 
the Euclidian norm is a natural choice in meshfree methods due to its rotational symmetry. Define 
the distance vectors Xi = Xi — xq Vi ~ 0, . . . ,m. The function value at each neighboring point 
u{xi) is given by a Taylor expansion as 

u{xi) = u{xq) + Vu(fo) • Xi + iV^u(xo) : {xi ■ xf) + . 

Here we use the Frobenius inner product A : B ~ Tli j '^ij^^ij- The error in the expansion is of 
order = 0{r^). A linear combination with coefficients (sq, . . . , Sm) equals 



Siu{xi) = u{xo) i'^Sij + Vu(fo) ■ ^ 



i=l / 
3\ 



This approximates the Laplacian, i.e. '^1^QSiu{xi) = Au{xq) + 0{r^), if exactness for constant, 
linear, and quadratic functions is satisfied 

mm m 

^ Sj = , ^ XiSi = , ^^{xi ■ xf) Si = 21 . (3) 

2=0 1 = 1 1=1 

Definition 1. A stencil (sq, ■ ■ • , Sm) to a set of points (a?0: ^^i; • • ■ i € B{xq, r) is called con- 
sistent (with the Laplace operator), if the constraints ([3]) are satisfied. 

The linear and quadratic constraints in ([3]) can be formulated as a linear system of equations 

V-s = b, (4) 

where V G R^x™ jg the Vandermonde matrix given by xi, . . . ,Xm, and s e M™ is the stencil 
vector. The number of constraints is fc = d{d+3} ^ rpj^^ constant constraint in ([3]) yields Sq = 
— ^i- Neumann boundary points can be treated in a similar manner. Approximating ^{xq) 

by X]i!lo ^i'^^i^i) leads to the constraints 

m m 

^ Si = , ^ XiS., = n . (5) 

1=0 i=l 
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For each point, a meshfree finite difference approximation consists of two steps: First, define which 
points are its neighbors. TypicaUy, more neighbors than constraints are chosen. Second, select a 
stencil. If ([3]) is underdetermined, a minimization problem is formulated to select a unique stencil. 
The results in Sect. [6] guarantee the existence of a solution to the above linear systems (|3|) and 

Definition 2. A consistent stencil {sq, . . . , Sm) is called minimal, if m < k. 

Minimal stencils are beneficial for the sparsity of the system matrix, resulting in a minimal memory 
consumption and a fast application of the matrix. In fact, the total number of neighboring points is 
proportional to the effort of multiplying the matrix with a vector. This promises a fast solution of 
system ^ with (semi-)iterative linear solvers, assumed the convergence properties are reasonable. 
More importantly, having minimal stencils results in an optimally sparse matrix graph. This 
property can be advantageous with respect to coarsening in algebraic multigrid methods. 

Definition 3. A consistent stencil {so,...,Sm) is called positive, if 
si, ■ ■ ■ ,Sm > 0. Due to ([3]) and ([5]), this implies for the central point sq < 0. 

Having positive stencils implies that the system matrix in ([2|) is an L-matrix (Def. [4|), which 
gives rise to an M-matrix structure (see Sect. SJ. The desirability of positive stencils has been 
pointed out by Demkowicz, Karafiat and Liszka [11]. Classical approaches do in general not yield 
positive stencils. An "optimal star selection" [TS] makes positive stencils likely, but they are not 
guaranteed. In Sect. [5] we present a strategy that approximates the Poisson equation ([1]) on a 
point cloud by minimal positive stencils. In Sect. [6] conditions on a point cloud are given which 
guarantee the existence of positive stencils. 

3. Non-Symmetric Matrices and Algebraic Multigrid 

Meshfree finite difference approaches, as described in Sect. [21 approximate the Poisson equation 
by non-symmetric matrices. Large linear systems with non-symmetric matrices have received con- 
siderable attention in Markov chain modeling. For instance, the computation of a page ranking in 
internet search engines leads to sparse non-symmetric matrices [26' . The arising matrices typically 
have an M-matrix structure. In the approximation of partial differential equations, non-symmetric 
matrices traditionally arise when convective terms are approximated. As described in [49], clas- 
sical AMG methods frequently are observed to converge for such non-symmetric matrices, if the 
advection does not dominate the problem. When approximating the Poisson equation, traditional 
approaches lead to symmetric matrices. In contrast, meshfree finite difference approaches approx- 
imate the Poisson equation by non-symmetric matrices. At first glance this sounds surprising, 
since the Laplace operator itself is selfadjoint. However, it can be proved and explained. 

The formal argument for non-symmetry is as follows. Consider two points Xi and Xj, each being 
a neighbor of the other, and a third point Xk which is a neighbor of Xi but not a neighbor of Xj. 
Since each stencil entry depends on all its neighbors, the point Xk influences the matrix entry Oy, 
but not the matrix entry Uji. 

The non-symmetry of the matrix can also be understood in that a self-adjoint operator is ap- 
proximated on a discrete set of points that is unstructured, and thus does not possess any local 
symmetry. Note that on unstructured geometries, classical finite element approaches [47] ap- 
proximate the Laplacian by a non-symmetric matrix as well. While in the finite element system 
Lil = M f, the stiffness matrix L alone is always symmetric, the actual approximation matrix 
M^^L is in general non-symmetric. 

The non-symmetry is inherent in the finite difference approximation, and to the author's knowl- 
edge, there is no straightforward way to symmetrize the linear system ([2]). This has interesting 
implications for the choice of efficient solvers for the arising linear systems. For instance, in the 
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realm of (preconditioned) Krylov subspace methods, the classical conjugate gradient method has 
to be replaced by more costly solvers, such as GMRES [H] or BiCGstab [5U] . 

Of fundamental interest are multigrid solvers for the arising systems ©. In most practical appli- 
cations with a sufficient number of particles, the largest part of the computational time is spent 
on the generation and solution of the linear systems in the pressure correction (and possibly in an 
implicit treatment of viscosity, which leads to quite similar system matrices). Hence it is of crucial 
interest to achieve the optimal computational effort 0{n), where n is the number of particles, that 
is promised by multigrid methods. 

In the context of particle methods, algebraic multigrid (AMG) approaches are most promising, 
since the application of geometric multigrid approaches on an unstructured, non-connected cloud 
of points is very difficult. While in the symmetric matrix case, various convergence results for 
AMG are well known ^ I48[ I49j , theoretical results for the non-symmetric matrix are much fewer 
in number. An analysis of AMG for non-symmetric matrices has been given by Notay |341 135j. A 
complementary analysis for methods of algebraic multilevel iteration (AMLI) type [U [2] has been 
given by Mense and Nabben [SUIISI]- AMLI methods are based on a block incomplete factorization 
of the matrix, partitioned in hierarchical form. Many AMG approaches can be interpreted as AMLI 
methods. 

A common assumption in all the above theoretical results is that the matrix be an M-matrix 
(see Sect. U]). The M-matrix structure implies certain properties that positive definiteness would 
imply for symmetric matrices. Since the negative Laplace operator in ([T]) is positive definite, the 
M-matrix structure is natural for approximations of it. In above references, proofs of convergence 
of two-grid methods for the M-matrix case are given; some results can be extended to the full 
multigrid case. In particular, it yields a discrete maximum principle, and guarantees the conver- 
gence of the Gauf5-Seidel iteration (see Sect. 14. 2p . In the above referenced theoretical results, the 
M-matrix structure is sufficient for the convergence of AMG and AMLI methods. The importance 
of the M-matrix structure for AMG is discussed in [34j . The analysis of the two-grid case is given 
in |35j . A key step in the generalization of theoretical AMG results to the non-symmetric matrix 
case is the suitable understanding of the size of an eigenvalue as its modulus in the complex case. 
In Sect. 14.31 recent theoretical results on the convergence of AMLI methods in the M-matrix case 
are outlined. 

Of course, the M-matrix property is not necessary for the convergence of multigrid approaches. 
In fact, convergence is frequently observed even if the M-matrix structure is violated. However, 
unless theoretical results for weaker assumptions are present, one cannot be sure. This would 
not be too problematic if a single Poisson equation had to be solved. However, if the Poisson 
solves are just substeps in a larger computation (as it is the case in the particle method), it is 
unsatisfactory to not have a guarantee of convergence. In fact, one may observe a rapid AMG 
convergence in the beginning of a long computation (when the particles are distributed nicely), 
but convergence may deteriorate as the particles become less evenly distributed. It is our hope 
that for meshfree finite difference methods, the M-matrix property can be a key instrument in 
successful and efficient solution approaches for the Poisson equation. In addition, it would be 
desirable to obtain an understanding of the connection between the performance of AMG and the 
corresponding particle geometry. 

4. Properties of M-Matrices 

In this section, we provide a brief overview of the definition of the M-matrix property, a practical 
sufficient condition that relates to meshfree finite difference approximations, and some well-known 
benefits of the M-matrix structure. 

Definition 4. A square matrix A = {aij)ij E M"^" is called Z-matrix, if Uij < ^ j ■ A 
Z-matrix is called L-matrix, if an > Vi. 
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We write A > if a.y > V«, j. The same notation applies to vectors. 
Definition 5. A regular matrix A is called inverse positive, if A^^ > 0. 
Definition 6. A Z-matrix is called M-matrix, if it is inverse positive. 

We use the M-matrix property, since it yields a sufficient condition for inverse positivity. 

4.1. A Sufficient Condition for the M-Matrix Structure. Conditions that imply the M- 
matrix property are required, since the inverse matrix is typically not directly available. Let the 
unknowns be labeled by an index set /. We consider square matrices A £ E^^^. 

Definition 7. The graph G(A) of a matrix A is defined by G{A) = {{i,j) e I x I : ^ 0}. The 
index i e / is called connected to j G /, if a chain i ~ i^^ii, . . . , i^-i, ik — j ^ I exists, such that 
) e G{A) Vi^ = l,...,fc. 

Remark 1. For a finite difference matrix, each index j G / corresponds to a point Xi. The index 
I G J being connected to j G / means that the point Xi connects (indirectly) to the point Xj via 
stencil entries. 

Definition 8. A finite difference matrix is called essentially irreducible if every point is connected 
to a Dirichlet boundary point. 

Remark 2. A finite difference matrix that is not essentially irreducible is singular, since the points 
that are not connected to a Dirichlet point form a singular submatrix. 

Definition 9. A matrix A G K^^^ is called essentially diagonally dominant if it is weakly diag- 
onally dominant (Vi G / : |aii| > X^fe^iti I'^ifeDi ^'^'^ every point i G / is connected to a point j £ I 
which satisfies the strict diagonal dominance relation \ajj\ > X^fc^^j l%fe|- 

Theorem 1. An L-matrix arising as a finite difference discretization of ^ is essentially diago- 
nally dominant if it is essentially irreducible. 

Proof. For an L-matrix the constant relation in ([3]) implies the weak diagonal dominance relation 
for every interior and Neumann point. Each row corresponding to a Dirichlet point satisfies the 
strict diagonal dominance relation. □ 

Theorem 2. An essentially diagonally dominant L-matrix is an M-matrix. 

Proof. The proof is given in [21, p. 153]. □ 

If problem ([1]) can be discretized by positive stencils and every point is connected to a Dirichlet 
point, then the resulting matrix is an M-matrix. 

4.2. Some Benefits of the M-Matrix Structure. The Poisson equation satisfies maximum 
principles. For instance, consider ([T]) with Dirichlet boundary conditions only. If / < and g < 0, 
then the solution satisfies u < [16]. A discretization by an M-matrix mimics this property in a 
discrete maximum principle. 

Theorem 3. Let A be an M-matrix. Then Ax < implies x < 0. Conversely, a Z-matrix 
satisfying Ax < ^ x < is an M-matrix. 

Proof. A is an M-matrix, thus A^^ > by definition. Let y — Ax. Then x — A^^y. The 
component-wise inequalities A^^ > and y < imply x < 0. The reverse statement is proved in 
[371 P- 29]. □ 

Theorem 4. If A is an M-matrix, and D its diagonal part, then p(I — D^^A) < 1, thus the Jacobi 
and the Gaufi-Seidel iteration converge. 
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Proof. The convergence of the Jacobi iteration is given in [3T]. The GauB-Seidel convergence 
foUows from the Stein- Rosenberg-Theorem [S2]- □ 

Further classical benefits of an M-matrix structure with respect to linear solvers can be found in 



4.3. Relevance of the M-Matrix Structure for Multilevel Methods. Recent theoretical 
convergence results for AMG in the non-symmetric M-matrix case have been given in [34|, I35j , 
and for AMLI type approaches in [30l [3T] . Here, we outline the key results for the latter type 
of approaches. The M-matrix structure is a sufficient condition for the convergence of additive 
Schwarz (AMLI) methods for non-symmetric matrices. Consider the system matrix be partitioned 
into blocks 

^ _ A-FF Afc 

\_AcF Acc\ 

where F denotes the fine grid unknowns, and C the coarse grid unknowns. If Aff is regular, then 
A can be factorized as 

A^\ " 
yAcFApp Ic 

where 

S = Acc - AcfApI,Afc (7) 
is the Schur complement. An approximation M ^ A is obtained by multiplying ([5]) with the 
approximations Aff ~ Aff and S ~ S. A simple example for S is to replace Aff by Aff in (H)- 
As outlined in [31], these block Schwarz factorization approaches are closely related to two-grid 
AMG methods. The AMLI method can be described as the stationary iteration with the iteration 
matrix 

T = I - M-^A . 

Definition 10. A splitting (M, N) oi A = M — N with M regular is called weak regular of the 
first type, if A/^^ > and M~^N > 0. 



Aff 
S 



If 




AppAFc 
Ic 



(6) 



The convergence of the AMLI method is ensured by the 

Theorem 5. // the splittings (^Aff, Aff — Aff^ o,nd (^S, S — are weak regular of the first 
type, then T > 0, and p(T) < 1. 



It is shown in |31j that if A is an M-matrix, then the assumptions in Thm. \5\ are satisfied for 
approximations Aff given by the Jacobi and the Gauss-Seidel methods and the incomplete LU 
factorization. However, the convergence proof fails in general, if A does not have an M-matrix 
structure. In [31], further methods have been presented (MAMLI, SMAMLI), that are based on 
other types of Schwarz factorizations. Under similar types of assumptions on the splittings, one 
has the convergence result p{Tsmamli) < p{Tmamli) < p{Tamli)- 



5. Least Squares vs. Linear Minimization Approach 

Classical approaches for meshfree derivative approximation are moving least squares methods, 
based on scattered data interpolation [53], and local approximation methods, based on generalized 
finite difference methods |28j . Their application to meshfree settings has been analyzed in [T^ [?7] . 

Around a central point aS'o, points inside a radius r are considered. A distance weight function 
w{5) is defined, which is small for 5 > r. For instance, one can consider w(5) = where a > 1. 
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Each point in the neighborhood Xi e B{xq, r) is assigned a weight Wi — w {\\xi — xqIU)- A unique 
stencil is defined via a quadratic minimization problem 

" 2 

min^^, s.t. y-s = 6. (8) 

Using = diag(wi), the solution of ([S]) is given by 

s^WV'^iVWV^y^ -h . (9) 

Least squares approaches do not yield minimal stencils, unless exactly k neighbors are considered. 
In general, they also do not yield positive stencils, as the following example shows. 

Example 1. Consider xq = (0,0) and 6 neighbors on the unit circle xi = (cos(-|(^i), sin(^iy9i)), 
where ((/?!, . . . , (^g) = (0, 1, 2, 3, 0.1, 0.2). Since all neighbors have the same distance from xq, the 
distance weight function does not play a role. Formula Q yields the non-positive least squares 
stencil s — (0.846, 1.005, 0.998, 1.003, 0.312, —0.164). However, the configuration admits a positive 
stencil, namely s = (1, 1, 1, 1, 0, 0). 

As outlined in Sects. [3land l4.3i we can be sure that multigrid approaches converge, if an M-matrix 
structure is generated. Hence, as a new approach, we enforce the positivity of stencils, i.e. we 
search for solutions in the polyhedron 

P = {s e M™ : s = &, s > 0} . (10) 

This is the feasibility problem of linear optimization 3l]. In Sect. [6] we provide conditions for the 
point geometry that guarantee the existence of solutions. If P is non-void and not degenerate, 
there are infinitely many feasible stencils. To single out a unique stencil we now formulate a linear 
(or minimization problem 

min V — , S.t. y • s = 6, s > , (11) 

^ Wi 

2—1 

where the weights Wi = w{\\xi — X0II2) are defined by an appropriately decaying (i.e. w(5) = 
o(|5|~^), see [45]) non-negative distance weight function w. Problem (jlip is a linear program (LP) 
in standard form. It is bounded, since we have imposed sign constraints and the weights Wi are 
all non-negative. 

Theorem 6. IJ the polyhedron (llOp is non-void, then the linear minimization approach (jlip yields 
minimal positive stencils. 



Proof. The sign constraints in (jlip ensure that the selected stencil is positive. The existence of 
a minimal solution is ensured by the fundamental theorem of linear programming |51j . If the LP 
(llip has a solution, then it also has a basic solution, in which at most k of the m stencil entries 
Si are different from zero. □ 



In contrast to least squares methods, the LP approach yields non-zero values only for a few selected 
points. This is in line with the work of Donoho |T3], who shows that a key advantage of 
minimization is sparsity. 

Remark 3. For regular grids, the linear minimization approach selects standard finite difference 
stencils. For instance, for a regular Cartesian grid, the classical 5-point stencil (2d) or 7-point 
stencil (3d), respectively, is obtained. In these cases, the basic solution is degenerate, i.e. some of 
the basis variables are zero. Also for the configuration given in Example[l] the linear minimization 
approach selects the classical 5-point Laplace stencil. 
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6. Existence of Positive Stencils and Matrix Connectivity 

The question of the existence of positive stencils, i.e. whether the polyhedron (fTO)) is non-void, is 
nontrivial. In [45] . the connection between the local point cloud geometry and the existence of a 
positive stencil has been presented. The key idea is to consider the dual linear program to (fTT|) by 
using Farkas' Lemma [51j . This yields a simple necessary criterion, as well as a simple sufficient 
criterion, as follows. 

Theorem 7 (Necessary Criterion for Positive Stencils). Let the considered central point w.l.o.g. lie 
in the origin. If a set of points X C around the origin admits a positive Laplace stencil, then 
they must not lie in one and the same half space ( with respect to an arbitrary hyperplane through 
the origin). 

Theorem 8 (Sufficient Criterion for Positive Stencils). Let the considered central point w.l.o.g. lie 

in the origin. Consider d G {2, 3}, and let (3 = \f2 — 1 in 2d, and (3 = y^i(3 — v^) in 3d. Assume 

the set of points X C M'* has the property that in any (i.e. v with \\v\\2 = 1 arbitrary) cone 

V ■ X > (l+/3^) ^|ja^||2; at least one neighboring point Xi £ X is contained. Then a positive 
Laplace .stencil exists. 

We call this criterion cone criterion. In 2d, the cone has a total opening angle 45°, and in 3d, 
the total opening angle is 33.7°. Using the cone criterion, the existence of positive stencils can be 
guaranteed by applying the £^ minimization (fTTj) to a large enough candidate set of neighboring 
points. As in [27], we define 

Definition 11. Let 51 C M'' be a domain and X — {xi, . . . , a?„} a point cloud. The mesh size h 
is defined as the minimal real number, such that fl C IJILi ^ (^^' l)^ where B {x,r) is the closed 
ball of radius r centered in x and Q is the closure of fl. 

Theorem 9 (Condition on Point Cloud to Guarantee Positive Stencils). Consider a point cloud 
of mesh size h, and let 7 be defined as in Thm. O // the radius of considered candidate points 
satisfies r > g^^(t^^2) f ' ^^^^^ f^^ every interior point which is sufficiently far from the boundary, a 
positive stencil exists. 

The proof, and a description of special conditions near the boundary, are given in [IS] . The ratio 
of candidate radius to maximum hole size radius equals 2.61 in 2d, and 3.45 in 3d. In practice, 
one generally observes that already much smaller smaller ratios yield positive stencils. 

Due to Thm. [1] and Thm. [21 the matrix composed of positive stencils is an M-matrix, if every 
interior and Neumann boundary point is connected to a Dirichlet boundary point. As discussed 
in [45j , the connectivity of every point to the boundary follows from Thm. [T] and the connectivity 
to a Dirichlet point can be ensured under comparably mild technical assumptions. Hence, if the 
implementation of the minimal positive stencil method ensures that Dirichlet boundary points are 
used in the stencils of nearby points, then the approach is guaranteed to generate M-matrices. 

7. Computational Effort to Generate the System Matrix 

In meshfree methods, the actual generation of the linear system that approximates the Poisson 
equation ([1]) carries a non-negligible computational cost. For each interior point, a minimization 
problem has to be solved. 

Classical least squares methods are based on the quadratic minimization ([5]) . The computation of 
the stencil by expression ^ requires first the set-up of the kx k matrix VWV^ , then the solution 
of the the linear system (V^VFF"^) -v = b, and lastly the computation of the product s = WV^ ■ v. 
This requires k{k + l)m + ^ floating point operations [121 p. 150]. 
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In contrast, the linear minimization approach (jlip does not admit an exphcit solution formula. 
Instead, a small linear program (LP) has to be solved. To our knowledge, there are no general 
results about efficient methods for such small LPs, especially any that would consider the special 
structure of the Vandermonde matrix. A numerical comparison of various methods has been 
presented in [43l p. 148]. Simplex methods [8] perform best for the arising LPs. A basis change 
corresponds to one stencil point replacing another. The theoretical worst case performance of 
simplex methods is not observed. Typical runs find the solution in about 1.5A: steps, resulting in 
a complexity of O(fc^m), which equals the asymptotic effort of least squares approaches. 

8. Numerical Results 

The two discretization approaches, least-squares vs. minimal positive stencils, have been com- 
pared in [45j in terms of accuracy of the finite difference approximation. For a sequence of point 
clouds, the global truncation errors with both approaches have been estimated numerically. The 
observation is that the approximation errors of the two approaches differ by at most 20%. An 
interesting aspect is that the measured order of convergence is generally second order, although 
the constraints ^ enforce only first order accuracy. This effect is also observed in other types of 
meshfree approaches. 

Since the two approaches do not differ significantly in the accuracy of approximation, a comparison 
of their computational effort when generating and solving the linear system is of interest. The 
presented linear minimization approach generates optimally sparse M-matrices, and we can hope 
that the solution of the arising linear systems is faster than for systems that are generated with 
classical least squares approaches. 

We consider two types of multigrid/multilevel approaches. In Sect. 18. ij we investigate the perfor- 
mance of SAMG ("Algebraic Multigrid Methods for Systems") [32] j by the Fraunhofer Institute 
for Algorithms and Scientific Computing. It is a library of AMG subroutines, based on the code 
RAMG05 |49l App. A], which is a successor of the pubhc domain code AMG1R5 [38]. SAMG is 
implemented as a preconditioner for a BiCGstab [50] iteration. It performs a full multigrid cycle 
with Gauss-Seidel relaxation, and sparse Gaussian elimination on the coarsest level. The numeri- 
cal tests presented here are performed using standard coarsening [321 PP- 473], though aggressive 
coarsening strategies [33 PP- 476] frequently turn out to improve the performance. In Sect. 18.21 
we consider an AMLI type method, implemented by C. Mcnse and R. Nabben, TU Berlin, 2005. 

8.1. Performance of SAMG. For a 3d Poisson test problem, the classical least squares approach 
is compared to the linear minimization approach in terms of the computational cost. We consider 
the geometry shown in Fig. [Tj It is motivated by the simulation of the flow of a viscous material 
in a geometry with a thin outlet. More detail is given in [43, p. 183]. The considered Poisson 
equation ([T]) has a smooth right hand side /. Neumann boundary conditions are prescribed ev- 
erywhere except for the thin outlet, at which Dirichlet boundary conditions are prescribed. For 
the considered geometry, the classical least squares approach is compared to the linear minimiza- 
tion approach in terms of the computational cost. For a sequence of point clouds with increasing 
resolution, the geometry is discretized. One approximation is generated using a weighted least 
squares approach ([8]), with about 40 neighbors for each interior point. Another approximation is 
generated from the £^ minimization approach (fTTjl . By construction, each interior point has only 
10 neighbors. Consequently, the £^ matrices possess about 4 times fewer non-zero entries than the 
LSQ matrices. 

Fig. [2| shows the CPU times for the setup of the system matrices. The solid curve represents the 
costs with the least squares approach, which is based on expression ([9]). The dashed curve shows 
the cost of the linear optimization approach, for which the simplex method is used. As predicted in 
Sect. [T] asymptotically the two approaches are equally expensive. In the considered example, the 
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Figure 1 . Geometry for the 3d Poisson test problem 




Figure 4. Memory consump- FIGURES. CPU times for solve 

tfon with SAMG with SAMG 



linear optimization approach is about 20% more expensive. Note that the implemented simplex 
algorithm is not optimized. Better linear programming methods, and/or smart heuristics may 
accelerate the matrix setup phase. 

Fig. [3] shows the CPU times for the solution of the arising system with a BiCGstab method [50] , 
preconditioned by ILU(O) [301 Chap. 10]. In comparison, Fig. O shows the CPU times for the 
solution of the same system with SAMG [42j. Three aspects can be observed: 

• The cost of the BiCGstab method grows superlinearly with the number of unknowns. In 
contrast, the SAMG cost is almost perfectly linear in the number of points. 
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• The SAMG method solves the arising Hnear systems significantly faster than the BiCGstab 
method. For low resolutions, SAMG is 2-3 times as fast as BiCGstab. For high resolutions, 
it is 4-5 times as fast. While the used BiCGstab method is certainly not fully optimized, 
the performance of SAMG is still compelling. 

• For both solution methods, the linear systems that arise from the £^ minimization approach 
are 3-4 times faster to solve than the systems generated by the LSQ approach. 

Fig. m shows the memory consumption in an SAMG solve. As one can expect, it is perfectly linear 
in the number of unknowns, and proportional to the number of non-zero matrix entries. 

The first conclusion we draw from this numerical experiment is that SAMG performs very well on 
both LSQ matrices and matrices. In particular, it is always observed to converge, even though 
the code is primarily designed for symmetric matrices, and even though most of the considered 
LSQ matrices are not M-matrices. SAMG is significantly faster than a classical BiCGstab method 
with ILU(O) preconditioning. In addition, the computational effort of SAMG is linear in the 
number of unknowns. 

The second conclusion is that the matrices obtained from the £^ minimization approach are solved 
significantly faster than classical LSQ matrices. The observed speed-up equals the increase in 
sparsity that the minimization approach yields. In another test (not shown here) , we store the 
non-basis neighbors, as obtained by the simplex method, as actual zero-entries in the matrix. The 
solution times with those matrices turn out to be within 10% of the solution times with the LSQ 
matrices. This implies that the M-matrix property itself does not increase the actual convergence 
rate of linear solvers. On the other hand, very sparse Poisson stencils are often times observed 
to be less robust than larger ones. In this sense, one can interpret the M-matrix property as an 
ingredient that allows the selection of a minimal stencil, without actually worsening convergence 
properties of linear solvers. 

8.2. Performance of an AMLI Type Method. We consider a 2d Poisson test problem on a 
circular domain = {a; e : ||x||2 < 1} as 

— Au = / in 
u = f on d^l 

where f{x, y) — sin(4a; + 0.1) + x cos(2j/ + 0.4). The analytical solution is u{x, y) — 16sin(4a; + 
0.1) -I- 4x cos(2j/-|-0.4). This problem is purposefully chosen to have a simple geometry, so that the 
focus lies on effects due to the unstructured point cloud. In [43l p. 173], this problem is used to 
analyze discretization errors on various point clouds. Here, we restrict to one specific point cloud, 
shown in Fig. [6l Boundary points are placed equidistantly on the unit circle. In the interior, 4000 
points are placed. For the given point cloud, we generate four types of discretization matrices: 
the £^ minimization matrix, and three LSQ matrices, with the closest 5/12/36 neighboring 
points considered in each stencil. For each matrix, the linear system is solved by various versions 
of AMLI type two-grid methods, all implemented by C. Mense and R. Nabben, TU Berlin, 2005. 

Two versions of coarsening are considered. In default coarsening, the first half of the unknowns 
is chosen as coarse scale variables, and the other ones are selected as fine scale variables. Since 
here the points are not sorted, this is equivalent to a random selection of coarse variables. In 
the Ruge-Stuhen coarsening [39], the coarsening of the matrix graph is performed preferably in 
the directions of large couplings, where the coupling strength (for an M-matrix) is defined as the 
magnitude of the corresponding off-diagonal entries. 

Four versions of AMLI type methods are considered, which are based on different splittings [30l[3T| . 
Specifically, besides the standard additive AMLI, a multiplicative (MAMLI), a reverse multiplica- 
tive (RMAMLI), and a symmetrized multiplicative (SMAMLI) splitting are considered. 
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Figure 6. Point cloud on Figure 7. Relative CPU times for various 

a circular domain types of AMLI approaches 



Fig. [7] shows the relativeQ CPU times for the various test cases. The left third corresponds to the 
£^ matrix, the middle third shows the LSQ approach with 12 neighbors for each point, and the 
right third represents the LSQ approach with 36 neighbors. The CPU times for the approach with 
the 5 nearest neighbors is not shown, since the AMLI methods failed to yield the correct solution. 
Needless to say, the obtained matrices are not M-matrices. However, also for the LSQ(12) and 
the LSQ(16) case, the Poisson matrices are in general not M-matrices. Nevertheless, the AMLI 
methods converge. In each case, the height of the dark part of the bar shows the setup time, and 
the light part shows the actual solution time. 

One can observe that the £^ matrices lead to the shortest solution times. However, unlike with the 
SAMG method, the speed-up factor with respect to the LSQ(12) case is not significant, although 
the £^ matrix is only half as dense. With respect to the LSQ(36) matrices, a noticeable speed-up is 
observed, although again not as significant as one would expect from the sparsity factors. Another 
observation is that for the considered example, there is no significant difference between the various 
splitting versions. Finally, for larger stencils, the Ruge-Stiiben coarsening is more expensive than 
the default coarsening. While this could be expected for the setup phase, it is surprising for the 
iteration phase. We cannot provide a satisfying explanation for this observation. 

9. Conclusions and Outlook 

With the pressure correction step in particle methods for incompressible fluid flows, we have 
presented an application that leads to large, sparse, non-symmetric matrices. These matrices that 
approximate the Poisson equation are constructed using a meshfree finite difference method. They 
have a different structure from traditionally considered non-symmetric matrices, such as Markov 
chain representation matrices. In particular, classical meshfree least squares approaches yield in 
general matrices that do not have an M-matrix structure. 

However, the M-matrix property is a key ingredient in recent convergence proofs of AMG and 
AMLI type methods for non-symmetric matrices. While it is not a necessary condition, conver- 
gence is not ensured if the M-matrix property is violated. In particle methods, Poisson equations 
are solved in every time step with ever changing geometries, and one single failure of convergence 
could spoil a whole long-term computation. Therefore, until multigrid convergence can be shown 
under weaker assumptions, the M-matrix structure is a property that should be guaranteed. The 
presented approach does so. Using an £^ minimization formulation, and linear programming meth- 
ods, optimally sparse M-matrices are generated. Conditions are provided that ensure the success 
of this approach. 



The absolute CPU times are not representative here and are thus omitted. 
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The performance of multigrid methods for meshfree finite difference approaches is numericahy 
investigated in various test cases. With SAMG, a sophisticated AMG method is apphed, and with 
AMLI, an approach is used in which convergence for M-matrices is proved. The main conclusion 
is that the minimization approach guarantees convergence and generates matrices that are 
significantly sparser than classical LSQ matrices, and this sparsity results in a speed-up in multigrid 
solvers. 

Currently, the setup of the system matrices with the £^ minimization approach is slightly slower 
than with LSQ approaches. However, more efficient linear programming methods may turn the 
tide towards the linear optimization approach. The efficient solution of the linear programs is a 
crucial point in the presented method, and worth a deeper analysis. Another interesting question is 
whether a weaker requirement than the M-matrix property already implies convergence of AMG 
and AMLI methods. Numerical tests presented in [33] indicate that the convergence of AMLI 
methods tends to be related to the inverse positivity of an LSQ matrix, rather than an M-matrix 
structure. In fact, LSQ matrices are frequently observed to be inverse positive, while not having 
an L-matrix structure. Alternative characterizations of inverse positive matrices |17] may help 
answer this question. 
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